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Abstract.  Discretization  of  the  Stokes  equations  produces  a  symmetric  indefinite  system  of  lin¬ 
ear  equations.  For  stable  discretizations,  a  variety  of  numerical  methods  have  been  proposed  that 
have  rates  of  convergence  independent  of  the  mesh  size  used  in  the  discretization.  In  this  paper,  we 
compare  the  performance  of  four  such  methods:  variants  of  the  Uzawa,  preconditioned  conjugate  gra¬ 
dient,  preconditioned  conjugate  residual,  and  multigrid  methods,  for  solving  several  two-dimensional 
model  problems.  The  results  indicate  that  where  it  is  applicable,  multigrid  with  smoothing  based 
on  incomplete  factorizaton  is  more  efficient  than  the  other  methods,  but  typically  by  no  more  than 
a  factor  of  two.  The  conjugate  residual  method  has  the  advantages  of  being  both  independent  of 
iteration  parameters  and  widely  applicable. 
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1.  Introduction.  Consider  the  system  of  partial  differential  equations 


(1) 


—A  u  +  Vp 
—  div  u  =  0 


u  =  0 
JnP  =  0 


/ 


on 

on  dCl, 


where  0  is  a  simply  connected  bounded  domain  in  Rd,  d  =  2  or  3.  This  system,  the 
Stokes  equations ,  is  a  fundamental  problem  arising  in  computational  fluid  dynamics, 
see  e.g.  [7,  12,  14,  17];  u  is  the  d- dimensional  velocity  vector  defined  on  0,  and  p 
represents  pressure. 

Discretization  of  (1)  by  finite  difference  or  finite  element  techniques  leads  to  a 
linear  system  of  equations  of  the  form 


(2) 


A  BT\U\=f 

B  —C  )  \  p  )  0 


where  A  is  a  set  of  uncoupled  discrete  Laplacian  operators  and  C  is  a  positive  semidef- 
inite  matrix.  We  consider  here  only  stable  discretizations,  i.e. ,  those  for  which  the 
condition  number  of  the  Schur  complement  matrix  BA~1BT  +  C  is  bounded  indepen¬ 
dently  of  the  mesh  size  used  in  the  discretization.  For  finite  element  discretizations 
with  C  =  0,  this  is  a  consequence  of  the  inf-sup  condition  and  upper  bound 


7  <  inf  sup 

<1  V 


(q,  div  v) 

Mi  Nlo  ’ 


|(g,  div  ^)| 
Milkllo 


where  7  and  T  are  independent  of  the  mesh  size.  Here,  |  •  |i  and  ||  •  ||o  denote  the 
ifCseminorm  and  Euclidean  norm,  respectively,  on  the  discrete  velocity  and  pressure 
spaces,  and  the  bounds  are  taken  over  all  v  and  q  in  the  appropriate  discrete  spaces; 
see  [7,  12,  14,  17]. 

In  recent  years,  a  variety  of  iterative  algorithms  have  been  devised  for  solving  the 
discrete  Stokes  equations.  In  this  paper,  we  compare  the  performance  of  four  such 
methods: 

1.  a  variant  of  the  Uzawa  method; 

2.  a  preconditioned  conjugate  gradient  (PCG)  method  applied  to  a  transformed 
version  of  (2); 

3.  a  preconditioned  conjugate  residual  (PCR)  method; 

4.  multigrid  (MG). 

The  Uzawa  method  is  the  first  among  these  to  have  been  devised  [2]  and  it  is  often 
advocated  as  an  efficient  solution  technique,  see  e.g.  [7,  12,  14].  The  convergence 
factor  associated  with  it  is  proportional  to  (k  —  1)/(k  +  1)  where  k  is  the  condition 
number  of  the  Schur  complement  BA_1BT  +  C  (see  §2.5).  The  conjugate  gradient 
method,  developed  by  Bramble  and  Pasciak  [5],  has  convergence  factor  proportional  to 
(\/k—  l)/(-v^Tl)  but  larger  cost  per  step  than  the  Uzawa  method.  The  preconditioned 
conjugate  residual  method  was  developed  by  Rusten  and  Winther  [24]  and  Silvester 
and  Wathen  [26,  31],  and  its  convergence  behavior  is  determined  by  properties  of  the 
indefinite  matrix.  For  multigrid,  we  consider  versions  derived  from  two  smoothing 
strategies:  a  variant  of  the  distributed  Gauss- Seidel  method  of  Brandt  and  Dinar  [6], 
and  the  technique  based  on  incomplete  factorization  developed  by  Wittum  [35]. 
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These  methods  all  have  the  property  that  for  appropriate  choice  of  precondition¬ 
ers  (or  for  multigrid,  smoothers),  their  convergence  rates  are  independent  of  the  mesh 
size  used  in  the  discretization.  The  actual  costs  of  using  them  depends  on  both  the 
convergence  rate  and  the  cost  per  iteration.  Our  goal  in  this  paper  is  to  compare  costs, 
in  operation  counts,  of  using  each  of  the  methods  to  solve  three  discrete  versions  of  (1). 
For  convergence  to  be  independent  of  mesh  size,  the  first  three  methods  ( Krylov  sub¬ 
space  methods )  require  a  preconditioning  operator  spectrally  equivalent  to  the  discrete 
Laplacian.  In  an  effort  to  unify  the  comparison  of  these  ideas  with  multigrid,  we  also 
implement  this  preconditioner  using  a  multigrid  method  for  the  associated  Poisson 
equation.  Our  main  conclusions  are  as  follows.  For  problems  where  it  is  applicable, 
one  version  of  multigrid,  using  incomplete  factorization,  requires  the  fewest  iterations 
and  operations,  but  it  is  only  marginally  faster,  i.e. ,  by  factors  of  approximately  1.5  to 
2,  than  the  Krylov  subspace  methods  and  the  distributed  Gauss-Seidel  method.  The 
Krylov  subspace  methods  are  more  widely  applicable  than  either  multigrid  method. 
Among  the  Krylov  subspace  methods,  the  conjugate  residual  method  is  slightly  slower 
than  the  conjugate  gradient  method,  and  in  some  cases,  the  Uzawa  method,  but  it 
has  the  advantage  of  not  requiring  any  parameter  estimates. 

An  outline  of  the  rest  of  the  paper  is  as  follows.  In  §2,  we  present  the  solution 
algorithms  and  give  an  overview  of  their  convergence  properties.  In  §3,  we  specify  four 
benchmark  problems  and  the  computational  costs  per  iteration  of  each  of  the  solution 
methods.  In  §4,  we  present  the  numerical  comparison. 

2.  Overview  of  methods.  In  this  section,  we  present  the  four  algorithms  un¬ 
der  consideration  and  outline  their  convergence  properties.  The  first  three  methods 
depend  on  a  preconditioning  operator  Qa  that  approximates  the  matrix  A  of  (2).  We 
assume  that  Qa  is  symmetric  positive  definite  (SPD)  and  that 


(3) 


Vi  < 


(v,  Av) 

(  ’  \  <  V2, 

{v,Qav) 


where  r] i  and  r] 2  are  independent  of  the  mesh  size  used  in  the  discretization.  In 
addition,  finite  element  discretizations  of  (1)  have  a  mass  matrix  M  associated  with 
the  pressure  discretization.1  The  preconditioner  will  also  include  a  SPD  approximation 
Qm  °f  M .  Discussions  of  computational  costs  will  be  made  in  terms  of  various  matrix 
operations  together  with  inner  products  and  “axpy’s,”  i.e.,  vector  operations  of  the 
form  y  ax  +  y. 

2.1.  The  inexact  Uzawa  method.  We  use  the  following  “inexact”  version  of 
the  Uzawa  algorithm  [11],  which  starts  with  uq  =  0  and  an  arbitrary  initial  guess  po: 


(4) 


for  i  =  0  until  convergence,  do 
Ui+ 1  =  ut  +  Q^if  ~  (Aui  +  BTpi)) 
pt+1  =  pt  +  a  QjJ(Bui+ 1  -  CpQ 

enddo 


Here,  is  a  scalar  parameter  that  must  be  determined  prior  to  the  iteration. 

In  the  “exact”  version  of  this  algorithm,  Qa  =  A  and  the  first  step  is  equivalent 
to  solving  the  linear  system  Aui+ 1  =  /  —  BTp{.  When  Qm  =  I,  the  exact  algorithm 

1  If  the  finite  element  solution  is  expressed  using  a  given  basis  {(j>t}  as  p  =  8,cj>t,  then  ||p||z2  = 

{8,  MS)1'2. 
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is  then  a  fixed  parameter  first  order  Richardson  iteration  applied  to  the  Schur  com¬ 
plement  system  ( BA~1BT  +  C)p  =  BA -1/;  Qm  is  a  preconditioner  for  this  iteration. 
The  inexact  Uzawa  algorithm  (4)  replaces  the  exact  computation  of  A_1(/  —  BTpi ) 
with  an  approximation. 

2.2.  A  preconditioned  conjugate  gradient  method.  Let  A  denote  the  co¬ 
efficient  matrix  of  (2).  Premultiplication  of  (2)  by  the  matrix 

Qa 1  o 

BQ~a 1  -I 

produces  the  equivalent  system 

(  Qa*  Qa*T  \(  u\  =  (  QaI  \ 

V  bQ~aa  ~  B  BQ~a1Bt  +  C  ){p  )  {  BQ-/J  )  ' 

Let  A 4  =  TA  denote  the  coefficient  matrix  of  this  system.  The  conjugate  gradient 
method  (CG)  developed  in  [5]  requires  that  the  bilinear  form 


(6) 


n 

qi 


(  v2 

\  ?2 


((A  -  Qa)v i,v2)  +  (qi,q2) 


define  an  inner  product.  Equivalently,  the  preconditioning  operator  Qa  must  satisfy 
(3)  with  r]i  >  1.  It  is  shown  in  [5]  that  Ai  is  SPD  with  respect  to  the  inner  product 
(6),  so  that  CG  in  this  inner  product  is  applicable.  The  matrix 


(7) 


Q 


I  0 
0  Qm 


is  also  SPD  with  respect  to  (6),  so  that  this  can  be  used  as  a  preconditioner. 
Let 


X0  = 


M0 

Po 


f  -  (Au0  +  BTp0) 
-( Bu0  -  Cpo ) 


denote  an  arbitrary  guess  for  the  solution  and  the  associated  residual.  An  implemen¬ 
tation  of  PCG  is  given  below.  Except  for  the  nonstandard  inner  product,  it  is  the 
standard  implementation,  as  given  for  example  in  [15,  p.  529].  It  is  more  efficient  then 
the  version  given  in  [5].  The  preconditioner  Qa  is  implicitly  incorporated  into  the 
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inner  product.  The  use  of  a  preconditioner  (7)  is  new. 

Ro  =  RRqi  Rq  =  Q  1Ro 
Po  =  Ro,  AiPo  =  P  APq 

a{0n)  =  [Ro,  Ro],  a{d)  =  [P0,MPo],  «o  =  a{Qn)/a{d) 

=  X0  +  op-Po 

R\  =  Ro  —  olqAPq  ,  i?i  =  i?o  —  aoA4Po,  R\  =  (7  ^R\ 

for  i  =  1  until  convergence,  do 

it  1  =  [Ri,  U  /3S  =  A-1  =  p£\ //3S 

Pi  =  Ri  +  Pi-lPi-l,  Ai  Pi  =  PAPi 

oh  =  dd.  o.«'‘>  =  [C„VIF,],  «,  =  «!*»/«!•» 

T  /' +  I  —  T /  T  ayT},  Ri-\-l  —  -R?  Oi^APi 

Ri+i  =  —  ay-AlPy  -Ri+i  =  t? 

enddo 

To  help  identify  operation  counts,  we  describe  the  computation  of  {ay}  and  {/3;} 
in  more  detail.  Letting 


Rt 


ri 


$i 


1 


R% 


ri 


> 


Ri 


we  have  =  [Rt,  Rt]  =  (rt,Art 


Si,  Si);  similarly,  if 


(8) 


Pi  = 


,  APi  — 


,  MPi 


then  af^  =  [Pi,  MPi]  =  (cy  AQ^Vi  —  Vi)  +  (p,  BQ~^Vi  —  ny).  Qp  is  referenced  only  in 
the  construction  of  Q~^v  in  (8),  so  that  only  the  action  of  the  inverse  of  Qa  is  required. 
Moreover,  although  the  vectors  Afi,  Aci  (for  ry)  and  AQ~^Vi  are  used,  the  first  two  of 
these  can  be  computed  using  an  AXPY.  Consequently,  only  one  matrix- vector  product 
by  A  is  needed. 

2.3.  The  preconditioned  conjugate  residual  method.  Since  A  is  symmet¬ 
ric,  variants  of  the  conjugate  residual  method  are  applicable.  Let  Xo  denote  the  initial 
guess  and  Rq  its  residual.  The  following  algorithm  implements  the  Orthomin  version 
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of  PCR  with  preconditioner  Q  [3]. 2 


Ro  =  Q-'Ro,  Po  =  Ro,  So  =  Q-MPo 

a{0n)  =  (Ro,AP0),  a{d)  =  (AP0,S0),  «o  =  a{Qn)/a{d) 

X\  =  Xo  +  aoPo,  R\  =  Rq  —  aoAPo,  R\  =  Rq  —  olqSq 
for  i  =  1  until  convergence,  do 

ft\  =  -(ARt,St.  i),  /3S  =  «a 

Pi  =  Ri  +  A-iP-1 ,  A  I)  =  5)  =  Q 

a|"}  =  (Ri,  APi),  af  =  (APi,  Si),  a =  of 

“t“  C^iPii  Ri-\- 1  —  Ri  ^ iAP{ ,  —  Ri 

enddo 

Any  symmetric  positive-definite  Q  could  be  used  as  a  preconditioner.  As  in  [26],  we 
use 


Qa  0 
0  Qm 


2.4.  Multigrid.  As  is  well  known,  multigrid  methods  combine  iterative  methods 
to  smooth  the  error  with  correction  derived  from  a  coarse  grid  computation.  We  use 
V-cycle  multigrid  for  “transformed  systems.”  Our  description  follows  [34,  35].  Cf. 
[22,  30]  for  other  multigrid  methods  derived  from  the  squared  system  associated  with 
(2). 

Let  —A p  denote  the  Laplace  operator  defined  on  the  pressure  space,  with  Neumann 
boundary  conditions  (see  [16]),  and  let  Ap  be  a  discrete  approximation  to  —  Ap  defined 
on  the  pressure  grid.  Consider  the  following  transformed  version  of  (2): 


(9) 


The  coefficient  matrix  in  (9)  is 

(10) 


A  = 


where  W  =  ABT  —  BT  Ap  and  G  =  BBT  -\-C Ap.  For  appropriate  discretizations  of  (1) 
(see  §3),  W  is  of  low  rank,  with  nonzero  entries  only  in  rows  corresponding  to  mesh 
points  next  to  90.  When  C  =  0,  G  can  also  be  viewed  as  discretization  of  —  Ap.  The 
splitting 


(11) 


A  =  s-n 


then  induces  a  stationary  iteration  applicable  to  (2), 


(12) 


Pk  +  1 


uk  |  ,  (  I  BT  \  x  (  f  -  ( Auk  +  BTpk ) 

Pk  )  l  0  -Ap  )  l  ~(Buk  -  Cpk) 


2  It  is  possible  for  this  version  of  PCR  to  break  down,  with  at  =  0.  The  ORTHODIR  version, 
which  uses  a  three-term  recurrence  to  generate  Pt,  is  guaranteed  not  to  break  down;  it  requires  two 
additional  AXPY’s.  Our  implementation  switches  from  the  ORTHOMIN  to  ORTHODIR  direction  update 
if  \at\  <  1 0 — 4 ,  as  described  in  [9].  In  the  experiments  discussed  in  §4,  this  switch  never  took  place. 
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This  is  used  as  the  smoother  for  the  multigrid  solver  for  (2).  Specific  choices  for  S  are 
given  in  §3.2. 

Let  Ru  denote  a  restriction  operator  mapping  velocity  vectors  in  the  fine  grid  (of 
width  h)  to  the  coarse  grid  (of  width  2 h),  let  Rp  similarly  denote  the  restriction  oper¬ 
ator  for  the  discrete  pressure  space,  and  let  Pu  and  Pp  denote  prolongation  operators 
from  the  coarse  spaces  to  the  fine  spaces.  (For  simplicity,  we  are  omitting  explicit 
mention  of  h  in  this  notation.)  One  step  of  V-cycle  multigrid  for  solving  (2),  starting 
with  initial  guess  u°,  p°,  is  as  follows. 


(u1,?1)  =  MG(u°,p°,f,g,k1,k2,h) 
if  h  <  ho,  then  %  Recursive  call 

Starting  with  u°,p° ,  perform  £q  smoothing  steps  (12),  producing  m1/3,^1/3 
rl/3  =  f_  (Au1/3  +  BTP'/3),  S !/3  =  —(Bu1/3  -  Cp1!3) 

r  c/3  =  Rur1/3,  sl/3  =  RpS1!3 

(u2J\pT)  =  MG(0,  0,  rc1/3,  sj/3,  ku  k2,  2 h) 


r2/ 3  =  u1!3  +  Pu 


2/3 


p2/3  _  pl/3  _|_  Ppp2J3 


Starting  with  n2/3,p2/3,  perform  k2  smoothing  steps  (12),  producing  u1  ,  p1 
else  %  Coarse  grid  solve  when  h  =  ho 


Solve 


A  Bt 
B  -C 


p 


f 

0 


directly 


endif 


We  also  use  V-cycle  multigrid  derived  from  the  discrete  Laplacian  as  a  preconditioner 
to  approximate  the  action  of  A~x  for  the  Krylov  subspace  methods;  this  is  defined 
analogously  and  we  omit  the  details.  For  all  multigrid  methods,  we  use  bilinear  inter¬ 
polation  to  define  Pu  and  Pp,  and  Ru  =  P(f ,  Rp  =  Pp  .  The  discrete  operators  at  each 
level  are  derived  from  the  discretization  on  the  associated  grid. 

2.5.  Convergence  properties.  We  briefly  outline  some  convergence  properties 
of  these  methods;  see  the  primary  references  for  derivations  of  bounds.  Each  of  the 
methods  generates  a  sequence  of  iterates  Ui  ~  u,  pi  ~  p  such  that,  if  e4-  is  a  represen¬ 
tation  of  the  error,  then  limj-^oodlej-H/HeoH)1/*  =  p  for  some  norm  ||  •  ||.  We  refer  to  p 
as  the  convergence  factor. 

We  are  assuming  that  the  discretization  and  choice  of  Qm  are  such  that 


(13) 


Ar  < 


(qABA-'B?  +  C)q) 

( q,QMq ) 


where  Ai  and  X2,  and  therefore,  k  =  X2/X2 ,  are  bounded  independently  of  the  mesh 
size  of  the  discretization.  This  is  the  case,  for  example,  when  Qm  is  a  suitable  ap¬ 
proximation  of  the  mass  matrix  in  finite  element  discretization  [29,  32].  Note  that  k 
is  the  spectral  condition  number  of  QfJ(BA~1BT  +  C). 

The  exact  Uzawa  algorithm  has  convergence  factor  p  (i  —  a  QfJ(BA~1BT  +  C)J 
[12].  This  is  smallest  for  the  choice  a  =  2/(Ai  +  A2),  in  which  case  it  has  the  value 
(k  —  1  )/(k+  1).  Thus,  the  convergence  factor  for  the  Uzawa  algorithm  is  independent 
of  the  mesh.  It  is  shown  in  [11]  that  the  performance  of  the  inexact  Uzawa  algorithm 
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is  close  to  that  of  the  exact  one  if  the  iterate  rq+i  satisfies 
(14)  ||/  -  BTpi  -  Aut+ 1 1|2  <  t\\ Bui  -  Cpi\ \Q-i  , 

where  r  is  independent  of  the  mesh  size. 

The  PCG  method  is  analyzed  in  [5,  Theorem  1],  where  it  is  shown  that  the  condi¬ 
tion  number  of  the  coefficient  matrix  Ai  of  (5)  is  bounded  by  a  constant  proportional 
to  k.  Thus,  standard  results  for  CG  [15]  imply  that  the  bound  on  the  convergence 
factor  for  this  method  is  proportional  to  (ffn  —  1  )/(\/k  +  1).  The  constant  of  pro¬ 
portionality  depends  on  how  close  r]  1  is  to  r] 2  in  (3),  i.e. ,  how  well  Qa  approximates 
A. 

The  PCR  method  is  analyzed  in  [24,  26].  The  analysis  shows  that  the  eigenvalues 
of  the  preconditioned  matrix  Q-1  A  are  contained  in  two  intervals  [—a,  —b]  U  [c,d], 
where  a,  b,  c,  are  d  are  positive  constants  that  are  independent  of  the  mesh  size.  The 
sizes  of  the  intervals  depend  on  k  and  the  accuracy  with  which  Qa  approximates  A.  It 
follows  from  the  convergence  analysis  of  CR  [9,  27]  that  the  convergence  factor  for  the 
preconditioned  algorithm  is  independent  of  the  mesh  size.  For  example,  it  is  shown 

[9]  that  if  d  —  c  =  a  —  6  >  0,  then  the  convergence  factor  is  bounded  by  2 

where  (3  =  (be) /(ad). 

It  is  shown  in  [36]  that  for  finite  difference  discretization  of  (1)  (see  §3.1),  two-grid 
variants  of  multigrid  are  convergent  with  convergence  rate  independent  of  the  mesh 
size.  The  analysis  applies  to  the  ILU  smoothing  of  §3.2,  although  it  requires  that  the 
prolongation  be  based  on  biquadratic  interpolation.  In  practice,  bilinear  interpolation 
has  been  observed  to  be  sufficient  [35].  Fourier  analysis  in  [6]  also  suggests  that 
MG/DGS  has  convergence  rate  independent  of  mesh  size. 

Remark  2.1.  Several  other  proposed  methods  share  properties  with  the  version  of 
PCG  under  consideration.  In  particular,  Verfiirth  [29]  has  shown  that  PCG  applied 
directly  to  the  Schur  complement  system  has  convergence  factor  proportional  to  pco] 
however,  this  method  requires  accurate  computation  of  the  action  of  A~x  at  each  CG 
step  [23].  Bank,  Welfert  and  Yserentant  [4]  present  a  method  making  use  of  Qa  ~  A 
with  convergence  rate  dependent  on  the  accuracy  of  this  approximation,  but  using  an 
additional  inner  iteration  on  the  pressure  space. 

3.  Solution  costs.  In  this  section,  we  outline  the  computational  costs  required 
to  solve  three  benchmark  problems  on  II  =  (0, 1)  X  (0, 1),  for  each  of  the  solution 
methods  of  §2. 

3.1.  Benchmark  problems.  We  use  four  discretizations  to  produce  test  prob¬ 
lems:  “marker  and  cell”  finite  differences,  and  three  mixed  finite  element  strategies. 

1.  Finite  differences  [19].  This  consists  of  the  usual  five-point  operator  for  each  of 
the  discrete  Laplacian  operators  of  (1),  together  with  centered  differences  for  the  first 
derivatives  Vp  and  div  u.  For  the  discretization  to  be  stable,  it  is  necessary  to  use 
staggered  grids  in  f l.  Figure  1  shows  such  grids  on  a  mesh  of  width  h  =  1/4.  In 
order  to  define  the  velocity  discretizations  at  grid  points  next  to  90,  certain  values 
outside  0  must  be  extrapolated;  for  example,  this  is  needed  to  approximate  d2ui/dy2 
for  points  “x”  next  to  the  bottom  of  90. 

2.  Linear/constant  finite  elements.  This  choice  consists  of  continuous  piecewise  linear 
velocities  on  a  mesh  of  width  h,  and  piecewise  constant  pressures  on  a  mesh  of  width 
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Fig.  1.  Staggered  grids  for  finite  difference  discretization. 
I~~  Cx)-|-  (x)— p-  j—  | 
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2/j.  The  discrete  pressures  are  not  required  to  be  continuous.  The  coarser  pressure 
grid  ensures  that  the  inf-sup  condition  holds  [17].  We  refer  to  this  as  the  Pi(h)Po(2h) 
discretization. 

3.  Piecewise  linear  finite  elements.  Here,  continuous  piecewise  linear  velocities  on 
a  mesh  of  width  h  are  paired  with  continuous  piecewise  linear  pressures  on  a  mesh 
of  width  2 h.  The  inf-sup  condition  is  also  satisfied.  We  call  this  the  Pi(h)Pi(2h) 
discretization. 

4.  Stabilized  piecewise  linear  finite  elements.  A  stable  discretization  using  piecewise 
linear  velocities  and  pressures  on  a  single  of  mesh  can  be  obtained  using  a  stabilization 
matrix  C  =  (3h2An ,  where  An  is  the  discrete  Laplace  operator  defined  on  the  pressure 
space,  subject  to  Neumann  boundary  conditions  [8].  This  technique  is  equivalent  to  the 
mini-element  discretization  [1]  after  elimination  of  the  internal  degrees  of  freedom.  We 
use  (3  =  .025,  as  recommended  in  [25].  We  refer  to  this  discretization  as  Pi(h)Pi(h). 
The  usual  hat  functions  are  used  as  the  bases  for  linear  velocities  and  pressures. 

The  coefficient  matrix  A  of  (2)  for  all  these  problems,  as  well  as  C ,  and 

BA~1Bt  +  C,  are  rank  deficient  by  one;  the  latter  three  matrices  share  a  constant 
null  vector.  As  a  result,  the  discrete  pressure  solutions  are  uniquely  defined  only  up  to 
a  constant.  In  exact  arithmetic,  the  solution  methods  under  consideration  correct  the 
initial  guess  with  quantities  orthogonal  to  the  null  space  of  A,  so  that  the  component 
of  the  null  space  in  the  computed  solution  is  the  same  as  in  the  initial  guess.  For  the 
analysis,  the  lower  bound  of  (13)  refers  to  the  smallest  nonzero  eigenvalue. 

Note  that  our  goal  in  considering  these  problems  is  to  compare  the  performance  of 
the  different  solution  strategies  on  a  variety  of  problems.  We  highlight  some  properties 
of  each  of  the  problems  as  follows: 

1.  finite  differences,  stable,  ^(pressure  unknowns)  ~  ^(velocity  grid  points); 

2.  finite  elements,  stable,  discontinuous  pressures,  ^(pressure  unknowns)  ~  ^ 
^(velocity  grid  points); 

3.  finite  elements,  stable,  continuous  pressures,  ^(pressure  unknowns)  ~  j  #(ve- 
locity  grid  points); 

4.  finite  elements,  requires  stabilization,  continuous  pressures,  ^(pressure  un¬ 
knowns)  ~  ^(velocity  grid  points). 

We  are  not  comparing  the  accuracy  achieved  by  the  discretizations,  and  remark  only 
that  the  three  finite  element  discretizations  display  the  same  asymptotic  convergence 
rates.  See  [17,  pp.  29,50]  for  comments  on  accuracy  of  finite  element  discretization, 
and  [21]  for  analysis  of  the  finite  difference  scheme. 

3.2.  Preconditioners  and  smoothers.  The  Uzawa,  PCR,  and  PCG  methods 
require  choices  of  Qa  and  Qm-  For  all  cases,  Qa  consists  of  one  step  of  V-cycle 


multigrid  derived  from  the  discrete  Laplacian.  To  ensure  that  Qa  is  symmetric,  the 
smoothing  is  based  on  damped  point- Jacobi  iteration  with  damping  parameter  uj  =  2/3 
[20].  For  the  three  finite  element  discretizations,  Qm  is  chosen  to  be  the  diagonal  of  the 
mass  matrix  iff,  see  [32].  (In  the  case  of  the  Pi(h)Po(2h)  discretization,  Qm  =  iff •) 
Although  there  is  no  mass  matrix  for  finite  differences,  a  natural  analogue  in  two 
dimensions  is  iff  =  h2I ,  and  this  is  used  for  Qm  with  finite  differences. 

We  consider  two  multigrid  smoothing  strategies.  The  first  is  a  variant  of  the 
distributed  Gauss-Seidel  (DGS)  iteration  introduced  by  Brandt  and  Dinar  [6].  The 
splitting  operator  of  (11)  is  given  by 


so  that  the  smoother  (12)  has  the  form 


uk+1  =  SA1(f  -  (Auk  +  BTpk)) 
pk+i  =  SQ1(-(B(uk  +  uk+1)  +  Cpk) 

uk+ 1  =  uk  +  Uk+1  +  BTPk+l 

Pk+1  —  Pk  Appk+1  ■ 


For  Sa ,  we  use  the  point  Gauss-Seidel  matrix  derived  from  red-black  ordering  of 
the  velocity  grid.  (That  is,  if  A  =  D  —  L  —  U  with  the  red-black  ordering,  then 
Sa  =  D  —  L.)  For  finite  differences,  Sq  =  (1/lj)T  where  T  is  the  tridiagonal  part 
of  G  and  uj  =  2/3;  that  is,  Sq  corresponds  to  a  damped  one-line  Jacobi  splitting. 
For  Pi(h)Pi(h)  finite  elements,  Sq  is  the  block  Jacobi  matrix  derived  from  a  two-line 
ordering  of  the  underlying  grid.  These  are  slightly  more  sophisticated  versions  of  the 
choice  Sq  =  diag(G)  used  in  [6].  We  refer  to  this  multigrid  method  as  MG/DGS. 

The  other  multigrid  smoother  is  the  incomplete  LU  factorization  (ILU)  presented 
by  Wittum  [35].  We  use  an  ILU  factorization  of  the  matrix  A  of  (10),  with  no  fill-in 
in  the  factors.  The  ordering  for  A  is  problem  dependent.  For  finite  differences,  it  is 
derived  from  an  uncoupled  red-black  ordering  of  the  underlying  grid.  That  is,  the  grid 
values  for  u\  were  listed  first,  in  red-black  ordering,  followed  by  those  for  U2,  and  then 
those  for  p.  (See  also  Remark  3.3  below.)  For  Pi(h)Pi(h)  finite  elements,  A  is  ordered 
according  to  an  uncoupled  lexicographic  ordering  of  the  grid  vectors.  We  denote  this 
method  by  MG/ILU. 

In  choosing  preconditioners  and  smoothers,  we  have  attempted  to  use  methods 
that  are  suitable  for  vector  and  parallel  computers.  Thus,  we  are  using  point  Jacobi 
smoothing  for  multigrid  preconditioning,  red-black  Gauss-Seidel  and  line  Jacobi  for 
the  DGS  iteration,  and  a  red-black  ordering  for  MG/ILU  applied  to  finite  differences. 
With  the  Pi(h)Pi(h)  discretization,  the  operator  G  in  the  DGS  method  is  a  19-point 
operator  that  has  block  Property  A  for  a  two-line  ordering  of  the  pressure  grid,  so 
that  the  two-line  Jacobi  splitting  can  be  implemented  efficiently  in  parallel.  The  ILU 
smoother  used  with  this  problem  is  not  efficient  on  parallel  computers.  Our  multigrid 
strategies  do  not  address  the  issue  of  idleness  of  parallel  processors  for  coarse  grid 
computations;  see  [10,  13]  for  discussions  of  this  point  for  the  discrete  Poisson  equation. 

Parameters  are  required  for  the  Uzawa,  PCG  and  multigrid  methods,  and  for  the 
multigrid  preconditioner.  These  are  as  follows: 
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Uzawa:  The  optimal  value  of  a  for  the  exact  Uzawa  method,  determined  empirically, 
is  used  for  the  inexact  version.  This  requires  computation  of  the  extreme 
eigenvalues  of  Q^(BA~1BT  +  C). 

PCG:  As  noted  in  §2.3,  the  preconditioner  must  be  scaled  so  that  771  >  1  in  (3). 
From  the  results  of  [5],  it  is  desirable  to  have  r]  1  close  to  1.  In  all  tests,  the 
scaling  is  chosen  so  that  1  <  r]i  <  1.02.  This  requires  computation  of  the 
smallest  eigenvalue  of  Q^A.3 

Multigrid:  For  the  coarse  mesh  size  ho  in  multigrid  computations,  we  chose  the  one 
of  ho  =  1/2  and  h0  =  1/4  that  produced  lower  iteration  counts.  This  turned 
out  to  be  ho  =  1/2  for  preconditioners  and  h0  =  1/4  for  solvers.  The  coarse 
grid  solution  is  obtained  using  Cholesky  factorization  for  the  preconditioners 
and  singular  value  decomposition  for  the  solvers. 

Remark  3.1.  For  the  Uzawa  method,  the  choice  of  Qa  does  not  guarantee  that  the 
condition  (14)  is  satisfied.  The  results  of  [11,  33]  as  well  as  those  of  §4  suggest  that 
with  multigrid  for  Qa ,  (14)  may  be  too  stringent. 

Remark  3.2.  The  effectiveness  of  the  multigrid  solvers  depends  on  the  fact  that  the 
commutator  W  in  (10)  is  zero  away  from  the  boundary  of  fh  This  is  true  for  the  finite 
difference  and  stabilized  Pi(h)Pi(h)  discretizations,  where  pressures  and  velocities  are 
defined  on  the  same  grid,  but  not  for  the  (stable)  Pi(h)Pi(2h)  discretization.  Our 
experiments  confirm  that  multigrid  is  ineffective  for  this  discretization,  and  we  do  not 
include  it  as  a  option.  See  [18,  p.  248]  for  a  discussion  of  this  issue.  For  the  Pi(h)Po(2h) 
discretization,  it  is  difficult  to  define  the  discrete  pressure  Poisson  operator  Ap,  and 
we  know  of  no  multigrid  implementation  for  this  problem. 

Remark  3.3.  For  MG/ILU  applied  to  the  finite  difference  discretization,  we  also 
tested  several  alternative  ordering  strategies,  including  an  uncoupled  lexicographic 
ordering  (i.e. ,  like  that  used  for  Pi(h)Pi(h)),  as  well  as  several  “coupled”  lexicographic 
orderings.  For  the  latter  strategies,  velocity  and  pressure  unknowns  are  not  separated 
from  one  another,  see  [28].  The  performances  of  MG/ILU  for  all  these  orderings  were 
very  close.  For  example,  for  h  =  1/32  as  in  Table  4  below,  the  smallest  average 
iteration  count  with  one  smoothing  step  was  10|  and  the  largest  was  ll|. 

3.3.  Iteration  costs.  We  identify  the  costs  per  iteration  of  each  of  the  methods 
by  first  specifying  the  “high  level”  operations  of  which  they  are  composed,  and  then 
determining  the  costs  of  each  of  these  operations.  High  level  operations  are  defined  to 
be  matrix- vector  products,  inner  products  (denoted  “(  ,  )”  in  the  tables  of  this  section), 
and  axpy’s.  Note  that  each  of  the  techniques  under  consideration  is  formulated  with 
essentially  the  same  set  of  these  operations;  consequently,  we  expect  operation  counts 
to  give  a  good  idea  of  their  comparative  performance. 

The  high  level  operations  are  shown  in  Table  1.  Matrix- vector  products  include 
operations  with  matrices  that  define  the  problem  or  method,  such  as  A  or  Ru ,  as  well 
as  preconditioning  and  smoothing  operators  such  as  Q^1  and  5//1.  The  latter  com¬ 
putations  are  themselves  built  from  other  matrix  operations,  and  some  of  these  are 
also  identified  in  the  table.  All  multigrid  entries  correspond  to  operations  performed 
on  one  grid  level.  For  multigrid  solvers,  the  smoothing  operations  are  presented  sep¬ 
arately;  these  operations  would  be  performed  £q  times  during  presmoothing  and 

3  In  the  experiments  described  in  §4,  these  were  computed  using  a  power  method  applied  to  Q~Q  A  — 
I;  five  to  ten  steps  were  needed  to  obtain  an  estimate  accurate  to  three  significant  digits. 
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Table  1 

High  level  operations  for  all  solution  algorithms. 


Matrix- Vector  Product 

AXPY 

(,  ) 

Uzawa 

1  A  1  Bt  1  Q~ax 
IB  1C  1QJ} 

1  ( nP ) 

1  (nu  +  np) 

PCG 

1  A  1  Bt  lQ^1 

2  B  1C  1QJ} 

4  (nu  +  np) 

2  (nu) 

3  (tiu  H-  ) 

PCR 

1  A  1  Bt  1QJ1 
IB  1C  1QJ} 

5  (nu  +  np) 

4  (nu  +  np) 

Multigrid 

Preconditioner 

(1  +  k\  +  ^2)  A  1  Ru 
(k\  +  ^2)  6/4 1  1  Pu 

Multigrid  Solver 

(Excluding 

smoother) 

1  A  1  Bt  1  Ru 

IB  1C  1  Rp 

1  Pu  1  Pp 

1  (nu  +  np) 

DGS  Smoother 

1  A  2  Bt  1  Ap 

IB  1C  U/1 

l^1 

ILU  Smoother 

1  A  2  Bt  1  Ap 
IB  1C  1  S-1 

times  during  postsmoothing.  The  lengths  of  the  vector  operations  are  listed  in  paren¬ 
theses.  We  are  assuming  that  one  inner  product  will  be  used  in  the  convergence  test, 
and  the  counts  in  the  table  include  this. 

The  costs  of  matrix- vector  products  are  estimated  to  be  the  number  of  nonzeros  in 
the  matrices  used.  This  is  roughly  one  half  the  number  of  “flops”  required,  and  it  is 
also  proportional  to  the  number  of  memory  references.  These  costs,  for  discretizations 
in  which  the  velocity  unknowns  come  from  an  n  X  n  grid,  are  shown  in  Table  2.  The 
costs  of  vector  operations  are  taken  to  be  the  length  of  the  vectors. 

Combining  the  data  of  Tables  1  and  2  gives  an  estimate  for  the  cost  per  iteration 
for  each  of  the  solution  methods  under  consideration.  These  numbers  are  all  propor¬ 
tional  to  n2,  and  we  present  in  Table  3  the  cost  factors  obtained  by  omitting  this 
factor,  rounded  to  the  nearest  integer.  For  the  multigrid  methods  (preconditioners 
and  solvers),  the  cost  of  one  full  multigrid  step  is  estimated  as  4/3  times  the  cost  of 
the  computations  on  the  finest  grid;  this  is  approximately  the  cost  of  full  recursive 
multigrid  in  two  dimensions. 

4.  Experimental  results.  We  now  present  the  results  of  numerical  experiments 
for  solving  (2).  All  experiments  were  performed  in  Matlab  on  a  Sparc-10  worksta¬ 
tion.  For  each  solution  algorithm,  we  solved  three  problems  derived  from  three  choices 
of  /  consisting  of  uniformly  distributed  random  numbers  in  [—1,1].  The  initial  guess 
in  all  cases  was  uq  =  0,  po  =  0.  The  stopping  criterion  was 

ll-Rilh/H-Rolh  <  10  6, 


n 


Table  2 

Costs  for  matrix-vector  products. 


Fin.  Diff. 

P1(h)P0(2h) 

P1(h)P1(2h) 

Pi(P)PiQ 

A 

10  ra2 

10ra2 

10  ra2 

B,  Bt 

4ra2 

4ra2 

8  ra2 

12ra2 

C 

0 

0 

0 

bn2 

Qm 

Ira2 

0.25ra2 

0.25ra2 

Ira2 

(Jacobi) 

2  ra2 

2  ra2 

2  ra2 

2  ra2 

(Gauss-Seidel) 

6ra2 

6ra2 

6ra2 

6ra2 

SB1 

3  ra2 

- 

- 

9  ra2 

Ap 

5  ra2 

- 

- 

bn2 

Pu  1  Pu 

6ra2 

4.5ra2 

4.5ra2 

4.5ra2 

Pp')  Pp 

3  ra2 

- 

- 

2.25ra2 

s -1 

19ra2 

41ra2 

Table  3 

Cost  factors. 

Uzawa  PCR  PCG  MG/DGS  MG/IC 

Finite 

Differences 

ki  =  k2  =  1 
ki  =  k2  =  2 

84  107  109  148  175 

116  139  141  244  297 

Pi(/0Po(2/0 

ki  =  k2  =  1 
k2  —  k2  —  2 

79  98  101 

111  130  133 

P1(/i)P1(2/i) 

ki  =  k2  =  1 
k2  —  k2  —  2 

86  104  111 

118  136  143 

Pi(h)Pi(h) 

k\  =  k2  =  1 
k2  =  k2  =  2 

101  124  134  247  333 

133  156  166  421  591 

Table  4 
Iterations. 


Uzawa  PCR  PCG  MG/DGS  MG/ILU 

Finite 

Differences 

k\  =  k2  =  1 

ki  =  k2  =  2 

36  41  30  24  12 

28  33  23  15  9 

P1(h)P0(2h) 

k\  =  k2  =  1 
k2  =  k2  =  2 

34  41  29  - 

26  34  23 

P1(h)P1(2h) 

k\  =  k2  =  1 
k2  —  k2  —  2 

89  57  38  - 

89  50  31 

Pi(h)Pi(fi) 

k\  =  k2  =  1 

ki  =  k2  =  2 

39  47  32  20  8 

38  40  25  10  7 

Table  5 

Estimates  of  convergence  factors. 


Uzawa  PCR  PCG  MG/DGS  MG/ILU 

Finite 

Differences 

k\  =  k2  =  1 

ki  =  k2  =  2 

.67  .70  .66  .62  .39 

.60  .64  .57  .50  .31 

P1(h)P0(2h) 

k\  =  k2  =  1 
k2  —  k2  —  2 

.69  .69  .70 

.58  .66  .55 

P1(h)P1(2h) 

k\  =  k2  =  1 
k2  —  k2  —  2 

.82  .79  .75 

.84  .78  .70 

Pi(h)Pi(fi) 

k\  =  k2  =  1 

ki  =  k2  =  2 

.70  .75  .68  .56  .24 

.70  .74  .62  .33  .21 

where 


We  found  that  performance  was  essentially  in  the  asymptotic  range  for  h  =  1/32,  and 
all  results  are  for  this  mesh  size. 

We  present  three  types  of  data:  iteration  counts,  estimates  for  convergence  factors, 
and  plots  of  residual  norms  as  functions  of  operation  counts.  The  iteration  counts  are 
averages  over  three  runs  of  the  number  of  steps  needed  to  satisfy  the  stopping  criterion; 
these  are  shown  in  Table  4.  The  estimates  for  asymptotic  convergence  factors  are  the 
averages  of  (H-Rs+il^/ll-RsIh)1^  over  all  steps  after  step  five;  here  Rk  represents  the 
average  of  the  fc’th  residual  norm  over  the  three  runs.  These  are  shown  in  Table  5. 
We  chose  step  five  rather  than  step  zero  because  performance  was  often  better  in  the 
first  few  steps  than  later,  when  the  asymptotic  behavior  is  seen.  Finally,  Figures  2  - 
5  plot  the  averages  of  the  residual  norms  against  operation  counts. 

We  make  the  following  observations  on  these  results. 

1.  Where  it  is  applicable,  multigrid  requires  the  smallest  number  of  iterations  and 
has  the  smallest  convergence  factors.  MG/ILU  is  superior  to  MG/DGS  in  these  mea¬ 
sures.  These  observations  agree  with  those  of  [35].  In  addition,  where  it  is  applicable, 
MG/ILU  requires  the  smallest  number  of  operations.  However,  multigrid  is  only  ef¬ 
fective  for  discretizations  where  velocities  and  pressures  are  defined  on  the  same  grid. 
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Fig.  2.  Operation  counts  for  finite  difference  discretization. 


Finite  Differences,  1  Smoothing  Step 


Fig.  3.  Operation  counts  for  Pi(h)Po('2h)  finite  element  discretization. 


2.  The  Krylov  subspace  methods  and  MG/DGS  are  roughly  equal  in  cost.  The  Krylov 
subspace  methods  are  more  widely  applicable  than  multigrid. 

3.  The  performances  of  all  these  methods  are  very  close.  In  terms  of  operation  counts, 
the  ratio  of  costs  of  the  most  expensive  and  least  expensive  method  is  no  worse  than 

2.3. 

4.  No  Krylov  subspace  method  is  clearly  superior  to  the  others.  PCG  exhibits  a 
somewhat  faster  convergence  rate  than  PCR,  and  the  Uzawa  algorithm  is  surprisingly 
competitive  with  the  other  two  methods.  This  appears  to  derive  from  the  dependence 
of  PCG  and  PCR  on  both  the  spectral  condition  number  k  from  (13)  and  the  accuracy 
of  the  preconditioner  Qa  as  an  approximation  to  A;  for  both  these  methods,  the 
iteration  counts  go  down  in  all  cases  when  the  number  of  smoothing  steps  in  Qa 
increases.  The  Uzawa  method  appears  to  be  less  sensitive  to  the  accuracy  of  Qa-  The 
values  of  k  for  the  three  problems  are: 

Finite  differences  4.14  P\(h)P\(2h)  22.71 

P1(h)P0(2h)  4.87  Pi(h)Pi(h)  9.91 

The  Uzawa  method  is  least  effective  for  the  Pi(h)Pi(2h)  discretization,  which  has  the 
largest  condition  number. 

5.  The  Uzawa  and  PCG  methods  depend  on  choices  of  iteration  parameters.  These  can 
be  estimated  relatively  inexpensively  (e.g.,  using  a  coarse  grid  for  the  Uzawa  method, 
and  a  few  steps  of  the  power  method  for  PCG),  but  this  increases  the  cost  of  these 
methods  and  makes  implementing  them  considerably  more  difficult.  In  contrast,  PCR 
is  independent  of  parameters  except  for  those  needed  for  the  multigrid  precondition¬ 
ing,  and  it  is  therefore  easier  to  implement.  Thus,  there  is  a  tradeoff  between  these 
methodologies:  PCR  converges  slightly  more  slowly  than  PCG  and,  often,  than  the 
Uzawa  method,  but  it  has  a  simpler  implementation. 

6.  For  each  of  the  solution  strategies  except  PCG,  it  is  less  expensive  to  use  one 
smoothing  step  than  two. 

Acknowledgements.  The  author  wishes  to  thank  David  Silvester  for  a  careful 
reading  of  a  preliminary  version  of  this  paper,  and  Andy  Wathen  for  some  helpful 
remarks. 


REFERENCES 

[1]  D.  Arnold,  F.  Brezzi,  AND  M.  Fortin,  A  stable  finite  element  for  the  Stokes  equations, 

Calcolo,  21  (1984),  pp.  337-344. 

[2]  K.  Arrow,  L.  Hurwicz,  AND  H.  Uzawa,  Studies  in  Nonlinear  Programming,  Stanford  Univer¬ 

sity  Press,  Stanford,  CA,  1958. 

[3]  S.  F.  Ashby,  T.  A.  Manteuffel,  AND  P.  E.  Saylor,  A  taxonomy  for  conjugate  gradient 

methods,  SIAM  J.  Numer.  Anal.,  27  (1990),  pp.  1542-1568. 

[4]  R.  E.  Bank,  B.  D.  Welfert,  AND  H.  Yserentant,  A  class  of  iterative  methods  for  solving 

saddle  point  problems,  Numer.  Math.,  56  (1990),  pp.  645-666. 

[5]  J.  H.  Bramble  AND  J.  E.  Pasciak,  A  preconditioning  technique  for  indefinite  systems  resulting 

from  mixed  approximations  of  elliptic  problems,  Math.  Comp.,  50  (1988),  pp.  1-17. 

[6]  A.  BRANDT  AND  N.  Dinar,  Multigrid  solutions  to  elliptic  flow  problems,  in  Numerical  Methods 

for  Partial  Differential  Equations,  S.  V.  Parter,  ed.,  Academic  Press,  New  York,  1979,  pp.  53- 
147. 

[7]  F.  Brezzi  AND  M.  Fortin,  Mixed  and  Hybrid  Finite  Element  Methods,  Springer- Verlag,  New 

York,  1991. 


16 


[8]  F.  Brezzi  AND  J.  PlTKARANTA,  On  the  stabilisation  of  finite  element  approximations  of  the 

Stokes  problem,  in  Efficient  Solutions  of  Elliptic  Systems,  W.  Hackbusch,  ed.,  Braunschweig, 
Vieweg,  1984,  pp.  11-19.  Notes  on  Numerical  Fluid  Mechanics,  Vol  10. 

[9]  R.  CHANDRA,  S.  C.  Eisenstat,  AND  M.  H.  Schultz,  The  modified  conjugate  residual  method 

for  partial  differential  equations,  in  Advances  in  Computer  Methods  for  Partial  Differential 
Equations  II,  R.  Vichnevetski,  ed.,  IMACS,  New  Brunswick,  1977,  pp.  13-19. 

[10]  N.  Decker,  Note  on  the  parallel  efficiency  of  the  Frederickson-McBryan  multigrid  algorithm, 

SIAM  J.  Sci.  Stat.  Comput.,  12  (1991),  pp.  208-220. 

[11]  H.  C.  Elman  AND  G.  H.  Golub,  Inexact  and  preconditioned  Uzawa  algorithms  for  saddle 

point  problems,  Tech.  Report  UMIACS-TR-93-41,  Institute  for  Advanced  Computer  Studies, 
University  of  Maryland,  1993.  To  appear  in  SIAM  J.  Numer.  Anal. 

[12]  M.  Fortin  AND  R.  Glowinski,  Augmented  Lagrangian  Methods:  Applications  to  the  Numerical 

Solution  of  Boundary-Value  Problems,  North-Holland,  New  York,  1983. 

[13]  P.  O.  FREDERICKSON  AND  O.  A.  McBryan,  Normalized  convergence  rates  for  the  PSMG 

method,  SIAM  J.  Sci.  Stat.  Comput.,  12  (1991),  pp.  221-229. 

[14]  R.  GLOWINSKI,  Numerical  Methods  for  Nonlinear  Variational  Problems,  Springer- Verlag,  New 

York,  1984. 

[15]  G.  H.  Golub  AND  C.  F.  Van  Loan,  Matrix  Computations,  The  Johns  Hopkins  University 

Press,  Baltimore,  second  ed.,  1989. 

[16]  P.  M.  Gresho  AND  R.  L.  Sani,  On  pressure  boundary  conditions  for  the  incompressible  Navier- 

Stokes  equations,  Int.  J.  Numer.  Meth.  Fluids,  7  (1987),  pp.  1111-1145. 

[17]  M.  GUNZBURGER,  Finite  Element  Methods  for  Viscous  Incompressible  Flows,  Academic  Press, 

San  Diego,  1989. 

[18]  W.  HACKBUSCH,  Multi-Grid  Methods  and  Applications,  Springer- Verlag,  Berlin,  1985. 

[19]  F.  H.  Harlow  AND  J.  E.  Welch,  Numerical  calculation  of  time-dependent  viscous  incompress¬ 

ible  flow  of  fluid  with  free  surface,  The  Physics  of  Fluids,  8  (1965),  pp.  2182-2189. 

[20]  S.  F.  McCormick,  ed.,  Multigrid  Methods,  SIAM,  Philadelphia,  1987. 

[21]  R.  A.  NlCOLAIDES,  Analysis  and  convergence  of  the  MAC  scheme  I,  SIAM  J.  Numer.  Anal.,  29 

(1992),  pp.  1579-1591. 

[22]  J.  PlTKARANTA  AND  T.  SAARINEN,  A  multigrid  version  of  a  simple  finite  element  method  for 

the  Stokes  problem,  Math.  Comp.,  45  (1985),  pp.  1-14. 

[23]  A.  Ramage  AND  A.  J.  WATHEN,  Iterative  Solution  Techniques  for  the  Navier-Stokes  Equations, 

Tech.  Report  93-01,  School  of  Mathematics  ,  University  of  Bristol,  1993.  To  appear  in  Int. 
J.  Numer.  Meth.  Fluids. 

[24]  T.  Rusten  AND  R.  WlNTHER,  A  preconditioned  iterative  method  for  saddle  point  problems, 

SIAM  J.  Matr.  Anal.  Appl.,  13  (1992),  pp.  887-904. 

[25]  D.  SILVESTER,  Optimal  low  order  finite  element  methods  for  incompressible  flow,  Comp.  Meths. 

Appl.  Mech.  Engrg.,  Ill  (1994),  pp.  357-368. 

[26]  D.  SILVESTER  AND  A.  Wathen,  Fast  Iterative  Solution  of  Stabilized  Stokes  Systems.  Part  II:  Us¬ 

ing  Block  Preconditioners,  Tech.  Report  218,  Mathematics  Dept.,  University  of  Manchester 
Institute  of  Science  and  Technology,  1992.  To  appear  in  SIAM  J.  Numer.  Anal.,  1994. 

[27]  D.  B.  SzYLD  AND  O.  B.  WlDLUND,  Variational  analysis  of  some  conjugate  gradient  methods, 

East- West  J.  of  Numer.  Math.,  1  (1991),  pp.  1-25. 

[28]  S.  P.  VANKA,  Block-implicit  multigrid  solution  of  Navier-Stokes  in  primitive  variables,  J.  Com¬ 

put.  Phys.,  65  (1986),  pp.  138-158. 

[29]  R.  VERFURTH,  A  combined  conjugate  gradient-multigrid  algorithm  for  the  numerical  solution  of 

the  Stokes  problem,  IMA  J.  Numer.  Anal.,  4  (1984),  pp.  441-455. 

[30]  - ,  A  multilevel  algorithm  for  mixed  problems,  SIAM  J.  Numer.  Anal.,  21  (1984),  pp.  264-271. 

[31]  A.  WATHEN  AND  D.  Silvester,  Fast  iterative  solution  of  stabilized  Stokes  systems.  Part  I:  Using 

simple  diagonal  preconditioners,  SIAM  J.  Numer.  Anal.,  30  (1993),  pp.  630-649. 

[32]  A.  J.  WATHEN,  Realistic  eigenvalue  bounds  for  the  Galerkin  mass  matrix,  IMA  J.  Numer.  Anal., 

7  (1987),  pp.  449-457. 

[33]  B.  D.  WELFERT,  Convergence  of  Inexact  Uzawa  Algorithms  for  Saddle  Point  Problems,  tech. 

report,  Mathematics  Department,  University  of  Arizona,  1993. 

[34]  P.  WESSELING,  An  Introduction  to  Multigrid  Methods,  John  Wiley  &  Sons,  New  York,  1992. 

[35]  G.  WlTTUM,  Multi-grid  methods  for  the  Stokes  and  Navier-Stokes  equations,  Numer.  Math.,  54 

(1989),  pp.  543-564. 

[36]  - ,  On  the  convergence  of  multi-grid  methods  with  transforming  smoothers,  Numer.  Math., 

57  (1990),  pp.  15-38. 


17 


